Response to a twist on systems with Z.p symmetry 
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We study response to a twist in the two-dimensional p-state clock model, which has the discrete 
Zp symmetry. The response is measured in terms of helicity modulus, which is usually defined with 
respect to an infinitesimal twist. However, we demonstrate that such a definition is inappropriate 
for the clock model. The helicity modulus must be defined with respect to a finite, quantized twist 
which matches the discrete Zp symmetry of the model. Recent numerical results, which casted doubt 
on the standard picture that two Berezinskii-Kosterlitz-Thouless transitions occur for p > 4, are 
rather a consequence of use of the inappropriate quantity. Numerical calculation of the appropriately 
defined helicity modulus indeed supports the standard picture. 
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Introduction. — Symmetries govern various phases of 
matter and transitions among them. There are two 
distinct classes of symmetries: continuous and discrete, 
with quite different consequences. For example, gapless 
Nambu-Goldstone modes are required in a spontaneous 
symmetry breaking phase, only if the symmetry is contin- 
uous. The cyclic group Zp symmetry is discrete for any 
value of p. Nevertheless, it approaches to the continuous 
U(l) symmetry in the limit of p — >■ oo. It is thus an in- 
teresting problem to investigate how the phase transition 
evolves as p is increased. 

Phase transitions in classical spin models, such as the 
classical XY model, in two dimensions with the contin- 
uous U(l) symmetry are an intriguing problem by itself. 
Berezinskii-Kosterlitz-Thouless (BKT) transition P,|2| is 
one of the most important concepts in statistical physics 
with a wide range of implications. The transition cannot 
be understood in terms of spontaneous symmetry break- 
ing. Yet, the response to a twist, which is quantified 
by the helicity modulus 3] , clearly distinguishes the two 
phases separated by the BKT transition: it is finite in 
the low-temperature (BKT) phase while it is zero in the 
high-temperature (disordered) phase. The helicity mod- 
ulus exhibits a universal jump at the BKT transition, 
which was confirmed numerically In the application 
to bosons in two dimensions, the helicity modulus is iden- 
tified with the superfiuid density. The universal jump of 
the helicity modulus implies the universal jump in su- 
perfiuid density which was remarkably observed in 
experiments 6]. 

Let us now discuss the case with the discrete Zp sym- 
metry. As an explicit example, we consider the p-state 
clock model on the two-dimensional L x L square lattice 
with periodic boundary conditions: 



Hp = — cos {Of' — 6pi — Ap^fr) . 



(1) 



(r,f') 



where r = (ri,r2) (ri_2 G Z, < ri_2 < L), {r,r') runs 
over all the nearest neighbor sites on the square lattice, 
and Of takes integral multiples of 27r/p. We have set 
the coupling constant to unity. The gauge field Ap^pi is 



usually set to zero but is introduced to impose a twist, 
as we will discuss later. 

According to the seminal renormalization group (RG) 
analysis , for p < 4, there is a single critical point sep- 
arating disordered and ordered phases, while the BKT 
phase appears for p > 4 in a finite range of temperature 
between the disordered and ordered phases. This has 
been also supported by later examinations ■ How- 

ever, several recent numerical studies pT| - [l4| challenge 
the standard picture. While details differ in each work, 
they suggest, in particular, the absence of the BKT tran- 
sition for p = 5. The most important evidence against 
the BKT transition picture is the lack of the universal 
jump of the helicity modulus at the transition. 

However, we will demonstrate that the apparent con- 
tradiction between the numerical results and the RG pic- 
ture is rather an artifact due to the helicity modulus de- 
fined inappropriately for the clock model. Once defined 
properly, the behavior of the helicity modulus is found to 
be perfectly consistent with the RG picture. 

Definition of helicity modulus for systems with dis- 
crete symmetry. — The helicity modulus is often de- 
fined with respect to an infinitesimal global twist. We 
first introduce the gauge field ^F-i-4,r — A/L for all 
the sites r, where x = (1,0) is the unit vector in the 
X direction, as a twist to the Hamiltonian ([1]). Then 
we define the helicity modulus by the second derivative 



d^Fp{A)/dA^^ 



A=0 



where Fp{A) is the free energy 



in the presence of the twist A introduced above. Being 
defined in terms of the derivative, it probes the response 
to an infinitesimal twist spread over the entire system. 

This definition of the helicity modulus is convenient 
for numerical calculations, since it can be written as an 
expectation value of a local physical quantity: 




cos 0p- 



r+x,r 




Sm 0f- 



r+x,r 



(2) 

where P — 1 jT is the inverse temperature (we normalize 
the Boltzmann constant to unity), and (jjp+x,? = Op+x — 
Op. We note that, for the infinitesimal global twist, the 
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resulting expression ([2]) is formally independent of p, and 
indeed the same as in the XY model (which corresponds 
to p — )■ oo). However, the actual value of course depends 
on p, since the expectation value is calculated for the p- 
state clock model in which the variable 9 is quantized to 
integral multiples of 27r/p. 

It was indeed this definition of the helicity modulus 
that was used in Refs. [ll|, [13] to study the p-state clock 
model. The most striking aspect of their results for p = 5 
is that, the helicity modulus thus obtained does not show 
the universal jump which is expected for a BKT transi- 
tion. In fact, the helicity modulus is nonvanishing even in 
the high-temperature disordered phase, in stark contrast 
to the expected behavior at a BKT transition. We also 
confirmed such behavior in our Monte Carlo simulations 
for p — 5. 

This seems to contradict with the standard RG pic- 
ture with two BKT transitions [7]. In the standard pic- 
ture, the helicity modulus should vanish in the disordered 
phase, while it diverges proportionally to the system size 
L in the ordered phase owing to the extra free energy of 
the induced domain wall. However, we find that Tp re- 
maining nonzero in the high-temperature phase is rather 
a consequence of the inappropriate definition of the he- 
licity modulus. Indeed, for any finite p, we shall demon- 
strate that the helicity modulus as defined in Eq. ^ is 
nonvanishing at arbitrary high temperature, where the 
system is certainly disordered, even in the thermody- 
namic limit. 

To see this, let us consider the high-temperature (char- 
acter) expansion of Eq. using the identity e'^'^°^'^ — 
T^V=-Qo-^k{P)e^^'^ for the modified Bessel function 1^. 
For the XY model, we find 
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,f,£Z]\r / {r,r') " 



(3) 



where kp pi is an integer defined for each link between 
the nearest neighbor pairs, and kpi p = —kp^p. The con- 
straint kp_p ~ means that the lattice divergence 
vanishes at every site. Thus, each term in the high- 
temperature expansion corresponds to a configuration of 
an integer- valued "flux field" {kf,fr} forming closed loops. 
Since /fe(/3) ~ Z?'*^', each term corresponds to the order 
J2{ff') \^r,r' \ of the high-tempcrature expansion. Unless 
there is a flux loop winding the entire system, which can 
occur only above the order L, we have '^frkpp^^ — 0. 
Thus, for the XY model in the thermodynamic limit 
L — > cx), the helicity modulus defined in Eq. ([2]) vanishes 
at any finite order in the high-temperature expansion ([3]). 

Now let us consider the p-state clock model. The only 
difference from the XY model comes from the fact that 
Ee=o,2^/p,...,2^(p-i) /p = P, if n is an integral muhiple 
of p. Thus, the high-temperature expansion of Eq. ^ for 
the p-state clock model is exactly the same as Eq. ([3]), 
except for the constraint at each site now replaced by 



^r,r' = mod p. This allows various flux configu- 
rations which were not permitted in the XY model. The 
lowest order in the high-temperature expansion among 
such configurations is given by fc^ — p and all other 
fiuxes being zero kp^p = 0. This leads to the nonvanish- 
ing contribution 



T 



P 



p-1 



(4) 



1)!2P 

in the lowest order of /?. This proves that the helicity 
modulus as defined in Eq. ([2]) does not vanish at arbitrary 
high temperature, for any finite p, even in the thermody- 
namic limit. This implies that Eq. ^ does not sharply 
distinguish the disordered and other phases. We note 
that, while this holds for an arbitrary large finite p, the 
contribution becomes rapidly smaller for larger p. This 
is the reason why the problem in the high-temperature 
phase has been noticed numerically only for small p, es- 
pecially p — 5. 

Thus we need a different quantity to describe the phase 
transition. Intuitively, the problem with the infinitesimal 
twist we have used could be understood as the introduc- 
tion of a mismatch between neighboring spins, whose di- 
rections are quantized in the p-state clock model. There- 
fore, let us define the helicity modulus as a response to 
the finite twist Ap^x,? = A localized on the horizontal 
links localized on a single column r G C = {(L — 1, r2)|0 < 
r-2 < L}: 

2{Fp{A) - Fp{0)) 



T — 



(5) 



Matching with the Zp symmetry of the model would re- 
quire A to be an integral multiple of 27r/p. We note that, 
in the original introduction of the helicity modulus [3], 
the twist was not restricted to be infinitesimal. Further- 
more, once inside the BKT phase, the helicity modulus 
thus defined should behave in the same way as the one 
with respect to an infinitesimal twist. 

In the high-temperature expansion, the difference be- 
tween Zp{A) and Zp{0) comes only from the factor 
exp [zA X^fec ^r.f+x] ■ Thus, if A is an integral multi- 
ple of 27r/p, the effect of the twist in the p-state clock 
model again disappears, even under the relaxed con- 
straint kp^pi = mod p, at any finite order of 
the high-temperature expansion in the thermodynamic 
limit. The helicity modulus thus defined indeed satisfies 
the physical expectation, vanishing exactly in the high- 
temperature disordered phase. 

Numerical calculations of the free energy difference due 
to a twist. — We have concluded that, in order to distin- 
guish different phases of the clock model in terms of its 
response to a twist, the twist angle must be an integral 
multiple of 27r/p. The change in the free energy due 
to a twist of a finite angle cannot be simply reduced to 
an expectation value of a local physical quantity, unlike 
in the case of the infinitesimal twist. Thus it is more 
difficult to calculate the response in Monte Carlo simula- 
tions. Nevertheless, it is possible to do so, as we demon- 
strate in the following. What we have to calculate is the 



3 




T 

FIG. 1. (Color online) (a),(b) Temperature dependence of Tp and Tp for p = 5 and p = 6. In both cases, Tp with respect to a 
infinitesimal twist does not reflect each phase appropriately. In comparison, for Tp with respect to the finite twist 2n /p, we can 
find three regions, blue, purple and beige regions, correspond to ordered, BKT and disordered phase, respectively. In addition, 
as L is increased, Tp approaches the universal value on the higher and lower transition points, with the universal relations ((B]) 
and (O shown respectively by yellow-green and aqua lines. Colored lines are for a guide to the eyes. (c),(d) Finite-size scaling 
of Ts near rgjfJf- and Tg^rp. All the data points collapse into universal curves in both cases. The estimated scaling parameters 
are (Tbkt, io) = (0.944, 0.8), (0.908, 1.9) for (c) and (d), respectively. 



ratio of the partition functions, Zp(A)/Zp(0), and this 
quantity can be calculated by boundary-flip Monte carlo 
method To overcome the critical slowing down of 

spin variables, we use Wolff cluster method [l3| at the low 
temperature. In our simulations, the number of Monte 
Carlo step are O(IO^) up to O(IO^) after equilibration. 

We first studied the helicity modulus ([5]) defined with 
respect to the finite quantized twist 27r/p for p < 4, 
where a single transition due to the spontaneous symme- 
try breaking is expected. We indeed find, as the system 
size is increased, Tp — in the high-temperature disor- 
dered phase and Tp oc L in the low-temperature ordered 
phase. The known critical temperatures for p = 2, 3, and 
4 are reproduced from the change in the behavior of Tp 
(data not shown in this Letter). 

Now let us come to the main question of p > 4, where 
the standard RG picture predicts two BKT transitions 
while some of the recent numerical studies found appar- 
ent contradictions with it. We emphasize that, while 
those recent papers do agree with the RG picture and 
identify BKT transitions for sufficiently large p, the he- 



licity modulus as defined in Eq. ([2]) they studied is non- 
zero at arbitrary high temperature. It is just that Eq. ([3]) 
is smaller for higher values of p and more difficult to be 
detected numerically. We still need to verify if the al- 
ternative definition, Eq. ([5]), indeed confirms the BKT 
transitions. We show the temperature dependence of the 
helicity modulus Tp, Eq. ([5]), with respect to the finite 
twist, for p = 5 and 6, in Fig. [IJa) and [Hb). In the 
low-temperature regime, colored blue, Tp increases pro- 
portionally to the system size L, while it decreases to- 
wards zero in the high-temperature regime, colored beige. 
These are precisely the expected behavior in the low- 
temperature ordered (spontaneously symmetry breaking) 
and the high-temperature disordered phases. In contrast, 
the helicity modulus Tp, defined as Eq. ^ with respect 
to the infinitesimal twist, does not exhibit physically ex- 
pected behaviors for a helicity modulus. While the dis- 
agreement of Eq. ^ with the standard RG picture in 
the high-temperature phase, especially for p — b, was 
discussed previously, we also emphasize that Eq. ([2]) does 
not detect the low-temperature ordered phase. The "fi- 
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nite twist" helicity modulus ([5]) we have introduced does 
distinguish the three phases, as is physicahy expected. 

Moreover, in the intermediate temperature regime 
(purple), Tp tends to be independent of L, which is in- 
deed the expected behavior in the BKT phase. This is 
rather clear for p = 6; the approximate transition tem- 
peratures estimated from Fig. [TJb) for p — 6, Tj 



0.90 and T^^j, 



high 
BKT 



timates T, 



high 



BKT 



~ 0.70, are in agreement with the es- 
0.90008(6) and T^^j, - 0.7014(11) ob- 
tained by a different method ^Ts'] . On the other hand, the 
existence of the BKT phase for p = 5 is not very clear 
in Fig. [IJa) , although it does not contradict the data. In 
addition, it is desirable to obtain more quantitative veri- 
fication also for p — 6. Therefore we further examine our 
numerical data with the BKT transition scenario. 

The well-known consequence of the BKT transition 
picture is the "universal jump" of the helicity modulus. 
In terms of RG, it is a consequence of the fact that the 
BKT transition occurs at a definite Gaussian coupling 
constant (also known as "Luttinger parameter" or "com- 
pactification radius" ) where the leading perturbation be- 
comes marginal. A similar argument can be also applied 
to the transition at the low-temperature side of the BKT 
phase f?\. For p > 4, we obtain the universal relations at 
the two transitions as 



lim 



-T, 



high 



BKT 



lim 



-r, 



n+0 



low 
BKT- 



(6) 
(7) 



In Fig. [Ija) andfTfb), we compare the temperature de- 
pendence of the "finite twist" helicity modulus (O, with 
the universal relations © and ([7]) shown respectively by 
yellow-green and aqua lines. We can find, as L is in- 
creased, Tp approaches the universal values on the higher 
and lower transition points. They do not, however, com- 
pletely converge owing to finite-size effects. 

Therefore, we have also performed a finite-size scal- 
ing analysis based on the RG theory at each transition. 
In the vicinity of the transition at higher temperature, 
-^BKT' anisotropy t/p is always irrelevant, so we can 

approximately neglect yp to obtain the standard BKT 
RG equations for the vortex fugacity y and the Gaussian 
coupling 0, 11k The scaling form of the helicity modulus 
was derived [l^ from the RG equations as 



x{T,L)^r'f{fd), 



(8) 



where x = 7r/3Tp — 2,6-/3 



/^BKT ^^'^ ^ ^ ln(i/Lo)- 
Here, /3bkt — 1/^bkt inverse transition tempera- 

ture, and Lq is a scaling parameter. They are the fitting 
parameters to the data. On the other hand, near the 
transition at the lower temperature, T^^r^, the vortex fu- 
gacity y is always irrelevant and thus may be neglected. 
The RG equations involving yp, in the vicinity of the 
lower-temperature transition, Tg^Ti similar to that 



for the higher-temperature transition at TgJlJ^ . As a con- 
sequence, for the lower-temperature transition at T'b°kT' 
the scaling equation is obtained from Eq. ([8]) by replacing 
X and 5 by x' = p"^ — 87r/?Tp and S' = /3bkt ~ respec- 
tively, where /3bkt = 1/^bkT' the scaling function 
/ by a different scaling function /'. 

In Fig. \V[c) and IHd), we show the finite-size scal- 
ing of T5 near each of the transitions, at T^kt ^^'^ 



T, 



low 
BKT- 



and TrTt- = 0.90514(9) estimated earher by a different 



We can find all the data points collapse into 
universal functions if we set two scaling parameters 
as (TBKT,io) = (0.944, 0.8), (0.908, 1.9), respectively. 
These are in good agreement with Tgjl^ = 0.95147(9) 

= 0.90514(9^ 
numerical method [20]. 

We have also confirmed, for p — 6, that Tg has the 
universal value at each BKT transition point. We esti- 
mate the transition temperatures Tg^J^ = 0.904(2) and 
^BKT = 0.700(2) by the finite-size scahng ^. These 
are m a very good agreement with T^'^^ = 0.90008(6) 
and Tgj^^ = 0.7014(11) determined earlier by a different 
numerical method |18| . 

Conclusions and discussions. — All our observations 
support, for p > 4, the existence of the BKT phase sepa- 
rated by two BKT transitions from the high-temperature 
disordered phase and from the low-temperature ordered 
phase, as predicted by RG Q. Much of the apparent con- 
tradiction with recent numerical results was rather due 
to the "helicity modulus" adopted in the numerical stud- 
ies, which was defined with respect to an infinitesimal 
twist. Such a quantity is inappropriate for the p-state 
clock model with the discrete Zp symmetry. We have 
established that, for this model, the helicity modulus in- 
stead has to be defined with respect to a finite quantized 
twist which matches the Zp symmetry of the system. 

Although we have restricted the analysis to the sim- 
ple p-state clock model in the present Letter, our discus- 
sion can be readily applied to more general models with 
Zp symmetry, including the XY model with a perturba- 
tion which breaks the U(l) symmetry down to Zp. The 
present result would also have implications in a wider 
range of problems with discrete symmetries or discretized 
degrees of freedom. 
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